#modifyId200712111134
#ADD output  path as parameter

#modifyId200712011838
#Renamed to UTR3ProEx
#3' utr extract

#modifyId20071201
#by zuofeng li
#modified for single sequence.

#this script is modified from GenbankToAnnotatedFasta.pl
#http://tiger.dbs.nus.edu.sg/microtar/files/GenbankToAnnotatedFasta.pl
#I reserved the origin function and add my own function
#zuofeng li
#Shanghai Center for  Bioinformation Technology
#2007.11.15

#!/usr/bin/perl
use strict;
use Bio::Seq;
use Bio::SeqIO;

{ # Start main block
	my $filename = shift || die("Paste a GenBank filename on the command line.\n\n");
	my $outpath = shift || die("Paste a output path without slash at the end on the command line.\n\n");
	my $i = 0;
	###
	my @fileList = ();	
	#####
	print "\nOpening source file $filename... ";
	my $seqIn = Bio::SeqIO->new(-format => 'GenBank', -file => "$filename");
	print "ok.\n";

	#for ($i = 0; my $input = $seqIn->next_seq; $i++) { #modifyId20071201
	my $input = $seqIn->next_seq;
		my $name = $input->display_id;
		my $sequence = lc($input->seq);
		my $description = $input->desc;
			$description =~ s/.$//;
		#my $id = "gi|". $input->primary_id ."|ref|". $name .".". $input->seq_version ."| ". $description;
		#modifyId20071201
		my $id = $input->primary_id;
		print "$name\n";

		my @cds = grep { $_->primary_tag eq 'CDS' } $input->get_SeqFeatures;
		foreach my $feature ( @cds ) {
			my $search = lc(($feature->seq)->seq);
			my $replace = uc($search);
			$sequence =~ s/$search/$replace/;
		}		
		my $outputfile = $outpath.$name.".fa";
		my $seqOut = Bio::SeqIO->new(-format => 'Fasta', -file => ">$outputfile");
		my $output = Bio::Seq->new(-id => $id, -seq => $sequence);
		$seqOut->write_seq($output);
#modifyId20071201
		#print $i;		
	#	@fileList[$i] = $name;
	#}
	#print "Processed ", $i, " input sequence(s).\n\n";
	
	####addd by zuofeng	
	print $name;
	&output3UTR($name, $outpath);
	print "Before go to next step, give a simple name to the UTR sequence!";
}

sub output3UTR()
{
	my ($name, $opath) = @_;
		unless(open(INFILE , "$name.fa"))
		{
			die("Open file input");
		}
		my $outfile = "3utr.fasta";
		
		my ($line);
		my ($seq);
		my $utrid;
		while($line=<INFILE>)
		{
			chomp($line);
			if($line=~/^>/){
				$line=~s/\>//;
				$utrid = $line;
			}
			else{
				$seq.= $line;
			}
		}
		
		my @field = split(//, $seq);
		
		my $init = 0;
		my $term = 0;
		my $k = 0;
		for($k =0; $k< $#field; $k++)
		{
			my $seg = $field[$k].$field[$k+1];
			#print $k." ".$seg."\n";
			if($seg=~/[A-Z][a-z]/){
				$term = $k+1;
			}
			
			if($seg=~/[a-z][A-Z]/){
				$init = $k +1+1;#the second "1" indicate the general coding system
			}
		}
		my $length = $k +1; 
		$utrid.= "|INIT|$init|Termination|$term|Length|$length";
		my $utrSeq = "";
		my $count =0;
		for($k=$term; $k<$length; $k++){
			$utrSeq.= $field[$k];
		}
		my $outurtpathfile = $opath.$outfile;
		my $seqOut = Bio::SeqIO->new(-format => 'Fasta', -file => ">$outurtpathfile");
		my $output = Bio::Seq->new(-id => $utrid, -seq => $utrSeq);
		$seqOut->write_seq($output);
		print "3'UTR sequence output to $outfile"."\n";
		close(OUTFILE);
		close(INFILE);
}